To run this vignette, we will use the Skin dataset from the article High-plex protein and whole transcriptome co-mapping at cellular resolution with spatial CITE-seq
To run this vignette, you must download the last version of Giotto Suite
remotes::install_github("drieslab/Giotto@suite_dev")
Load Giotto
library(Giotto)
## Giotto Suite 3.3.1
instrs <- createGiottoInstructions(save_plot = FALSE,
show_plot = FALSE)
Create Giotto object using RNA and Protein expression, as well as spatial positions.
my_giotto_object <- createGiottoObject(
expression = list(rna = list(raw = "data/expression_rna.csv"),
protein = list(raw = "data/expression_protein.csv")),
expression_feat = list("rna", "protein"),
spatial_locs = "data/positions.csv",
instructions = instrs)
my_giotto_image <- createGiottoImage(gobject = my_giotto_object,
do_manual_adj = TRUE,
scale_factor = 0.0625,
mg_object = "data/skin.jpg",
negative_y = TRUE)
## do_manual_adj == TRUE
## Boundaries will be adjusted by given values.
my_giotto_object <- addGiottoImage(gobject = my_giotto_object,
images = list(my_giotto_image),
spat_loc_name = "raw")
Visualize image
spatPlot2D(my_giotto_object,
show_image = TRUE,
point_size = 2)
my_giotto_object <- filterGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_threshold = 1,
feat_det_in_min_cells = 1,
min_det_feats_per_cell = 1)
## completed 1: preparation
## completed 2: subset expression data
## completed 3: subset spatial locations
## completed 4: subset cell (spatial units) and feature IDs
## completed 5: subset cell metadata
## completed 6: subset feature metadata
## completed 7: subset spatial network(s)
## completed 8: subsetted dimension reductions
## completed 9: subsetted nearest network(s)
## completed 10: subsetted spatial enrichment results
##
## Feature type: rna
## Number of cells removed: 0 out of 1691
## Number of feats removed: 363 out of 15486
my_giotto_object <- normalizeGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
norm_methods = "standard",
scalefactor = 10000,
verbose = TRUE)
##
## first scale feats and then cells
my_giotto_object <- addStatistics(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna")
my_giotto_object <- calculateHVF(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_values = "normalized",
method = "cov_groups",
nr_expression_groups = 20,
zscore_threshold = 1.5)
## return_plot = TRUE and return_gobject = TRUE
##
## plot will not be returned to object, but can still be saved with save_plot = TRUE or manually
my_giotto_object <- runPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_values = "normalized",
reduction = "cells",
feats_to_use = "hvf",
name = "rna.pca")
## "hvf" was found in the feats metadata information and will be used to select
## highly variable features
## class of selected matrix: dgCMatrix
## [1] "finished runPCA_factominer, method == factominer"
my_giotto_object <- runUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_values = "normalized",
reduction = "cells",
dimensions_to_use = 1:10,
dim_reduction_name = "rna.pca")
my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
type = "kNN",
dim_reduction_name = "rna.pca",
name = "rna_kNN.pca",
dimensions_to_use = 1:10,
k = 20)
my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
nn_network_to_use = "kNN",
network_name = "rna_kNN.pca",
name = "leiden_clus")
plotPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
dim_reduction_name = "rna.pca",
cell_color = 'leiden_clus',
title = "RNA PCA")
plotUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
cell_color = 'leiden_clus',
point_size = 2,
title = "RNA Uniform Manifold Approximation & Projection (UMAP)",
axis_title = 12,
axis_text = 10)
spatPlot2D(my_giotto_object,
show_image = TRUE,
point_size = 1.5,
cell_color = "leiden_clus",
title = "RNA Leiden clustering")
my_giotto_object <- filterGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_threshold = 1,
feat_det_in_min_cells = 1,
min_det_feats_per_cell = 1)
## completed 1: preparation
## completed 2: subset expression data
## completed 3: subset spatial locations
## completed 4: subset cell (spatial units) and feature IDs
## completed 5: subset cell metadata
## completed 6: subset feature metadata
## completed 7: subset spatial network(s)
## completed 8: subsetted dimension reductions
## completed 9: subsetted nearest network(s)
## completed 10: subsetted spatial enrichment results
##
## Feature type: protein
## Number of cells removed: 0 out of 1691
## Number of feats removed: 0 out of 283
my_giotto_object <- normalizeGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
scalefactor = 10000,
verbose = T)
##
## first scale feats and then cells
my_giotto_object <- addStatistics(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_values = "normalized")
my_giotto_object <- runPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_values = "normalized",
scale_unit = T,
center = F,
method = "factominer")
## "hvf" was not found in the gene metadata information, all genes will be used
## class of selected matrix: dgCMatrix
## [1] "finished runPCA_factominer, method == factominer"
my_giotto_object <- runUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_values = "normalized",
dimensions_to_use = 1:10)
my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
type = "kNN",
name = "protein_kNN.pca",
dimensions_to_use = 1:10,
k = 20)
my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
nn_network_to_use = "kNN",
network_name = "protein_kNN.pca",
name = "leiden_clus")
plotPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
dim_reduction_name = "protein.pca",
cell_color = 'leiden_clus',
title = "Protein PCA")
plotUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
cell_color = 'leiden_clus',
dim_reduction_name = "protein.umap",
point_size = 2,
title = "Protein Uniform Manifold Approximation & Projection (UMAP)",
axis_title = 12,
axis_text = 10)
spatPlot2D(my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
cell_color = "leiden_clus",
point_size = 1.5,
show_image = TRUE,
title = "Protein Leiden clustering")
saveGiotto(my_giotto_object, "multiomics_Giotto_object")
Download the file from https://drive.google.com/file/d/1GyhKWl14MxEJCw_2-e8Bge3d3zY7XLkc/view?usp=share_link
my_giotto_object <- loadGiotto("multiomics_Giotto_object")
my_giotto_object <- runWNN(my_giotto_object,
modality_1 = "rna",
modality_2 = "protein",
pca_name_modality_1 = "rna.pca",
pca_name_modality_2 = "protein.pca",
k = 20,
verbose = TRUE)
## The NN network name was not specified, default to the first: "rna_kNN.pca"
## The NN network name was not specified, default to the first: "protein_kNN.pca"
## [1] "Calculating rna - rna distance"
## [1] "Calculating protein - protein distance"
## [1] "Calculating low dimensional cell-cell distances for rna"
## [1] "Calculating low dimensional cell-cell distances for protein"
## [1] "Calculating within-modality prediction"
## [1] "Calculating cross-modality prediction"
## [1] "Calculating Jaccard similarities"
## [1] "Calculating kernel bandwidths"
## [1] "Calculating modality weights"
## [1] "Calculating WNN"
## [1] "Calculating WNN normalization"
## [1] "Calculating WNN graph"
## [1] "Saving WNN results"
my_giotto_object <- runIntegratedUMAP(my_giotto_object,
modality1 = "rna",
modality2 = "protein",
spread = 5,
min_dist = 2)
## [1] "Calculating integrated Nearest Neighbors"
## [1] "Creating igraph"
## [1] "Calculating integrated UMAP"
my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
nn_network_to_use = "kNN",
network_name = "integrated_kNN",
name = "integrated_leiden_clus",
resolution = 0.4)
plotUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
cell_color = 'integrated_leiden_clus',
dim_reduction_name = "integrated.umap",
point_size = 1.5,
title = "Integrated UMAP using Integrated Leiden clusters",
axis_title = 12,
axis_text = 10)
spatPlot2D(my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
cell_color = "integrated_leiden_clus",
point_size = 1.5,
show_image = TRUE,
title = "Integrated Leiden clustering")
my_giotto_image <- createGiottoImage(gobject = my_giotto_object,
do_manual_adj = TRUE,
scale_factor = 0.0625,
mg_object = "data/skin.jpg",
negative_y = TRUE)
my_giotto_object <- addGiottoImage(gobject = my_giotto_object,
images = list(my_giotto_image),
spat_loc_name = "raw")
The interactive selection tool will allow you to run a Shiny app and plot polygons over the tissue.
Note that you can use the blue bars to Zoom in/out at areas of interest over the tissue
You can choose multiple polygon names, one for each drawn polygon.
When you have finished plotting polygons, click on the Done button
my_polygon_coordinates will store a data.table object with the polygon names and xy coordinates for downstream analysis.
my_spatPlot <- spatPlot2D(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
show_image = TRUE,
cell_color = 'leiden_clus',
point_size = 2,
point_alpha = 0.75)
my_spatPlot
my_polygon_coordinates <- plotInteractivePolygons(my_spatPlot, color = "black")
my_giotto_polygons <- createGiottoPolygonsFromDfr(my_polygon_coordinates,
name = 'selections')
## Selecting col "name" as poly_ID column
## Selecting cols "x" and "y" as x and y respectively
my_giotto_object <- addGiottoPolygons(gobject = my_giotto_object,
gpolygons = list(my_giotto_polygons))
addPolygonCells will modify the cell_metadata slot to indicate if a certain cell is located within a polygon area (polygon1,2,3..) or not (no_polygon)
my_giotto_object <- addPolygonCells(my_giotto_object,
polygon_name = 'selections')
##
## These column names were already used: nr_feats perc_feats total_expr
## leiden_clus integrated_leiden_clus
## and will be overwritten
Verify that cell metadata has been modified
head(pDataDT(my_giotto_object))
## cell_ID nr_feats perc_feats total_expr leiden_clus integrated_leiden_clus
## 1: 10x14 70 0.4628711 474.8438 2 7
## 2: 10x15 237 1.5671494 1228.7768 4 7
## 3: 10x16 92 0.6083449 599.3736 2 5
## 4: 10x17 87 0.5752827 581.0906 5 4
## 5: 10x18 132 0.8728427 807.0911 7 1
## 6: 10x19 336 2.2217814 1601.3652 3 1
## selections
## 1: polygon 2
## 2: polygon 2
## 3: polygon 2
## 4: polygon 2
## 5: no_polygon
## 6: no_polygon
my_polygons_cells is a terra::spatVector object that contains the information of each cell within the polygon areas
my_polygons_cells <- getCellsFromPolygon(my_giotto_object,
polygon_name = 'selections')
my_polygons_cells
## class : SpatVector
## geometry : points
## dimensions : 465, 4 (geometries, attributes)
## extent : 1, 50, -43, -1 (xmin, xmax, ymin, ymax)
## coord. ref. :
## names : sdimx sdimy cell_ID poly_ID
## type : <int> <int> <chr> <chr>
## values : 10 -14 10x14 polygon 2
## 10 -15 10x15 polygon 2
## 10 -16 10x16 polygon 2
Find marker features for each polygon
scran_results <- findMarkers_one_vs_all(my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
method = "scran",
expression_values = "normalized",
cluster_column = "selections",
min_feats = 10)
## using 'Scran' to detect marker feats. If used in published research, please cite:
## Lun ATL, McCarthy DJ, Marioni JC (2016).
## 'A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor.'
## F1000Res., 5, 2122. doi: 10.12688/f1000research.9501.2.
##
## start with cluster no_polygon
##
## start with cluster polygon 1
##
## start with cluster polygon 2
top_genes <- scran_results[, head(.SD, 2), by = 'cluster']$feats
comparePolygonExpression will create a heatmap to observe the differential expression between polygon areas
comparePolygonExpression(my_giotto_object,
selected_feats = "top_genes")
## using 'Scran' to detect marker feats. If used in published research, please cite:
## Lun ATL, McCarthy DJ, Marioni JC (2016).
## 'A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor.'
## F1000Res., 5, 2122. doi: 10.12688/f1000research.9501.2.
##
## start with cluster no_polygon
##
## start with cluster polygon 1
##
## start with cluster polygon 2
You can isolate each region to plot cell metadata such as the number of features per cell.
spatPlot2D(my_giotto_object,
cell_color = 'nr_feats',
group_by = 'selections',
color_as_factor = FALSE,
point_size = 2,
coord_fix_ratio = 1,
cow_n_col = 2,
show_legend = TRUE)
Or plot gene expression per cell by selected region (polygon).
spatFeatPlot2D(my_giotto_object,
expression_values = 'scaled',
feats = c("APOC1", "TMEM132D"),
group_by = 'selections',
point_size = 2)
You can compare the percent of cell types by selected region (polygon).
compareCellAbundance(my_giotto_object,
spat_unit = "cell",
feat_type = "rna")
You can isolate each region and plot cell types
## cell_ID selections
## 1: 10x14 polygon 2
## 2: 10x15 polygon 2
## 3: 10x16 polygon 2
## 4: 10x17 polygon 2
## 5: 10x18 no_polygon
## ---
## 1687: 9x46 no_polygon
## 1688: 9x47 no_polygon
## 1689: 9x48 no_polygon
## 1690: 9x49 no_polygon
## 1691: 9x50 no_polygon
spatPlot2D(my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
cell_color = 'integrated_leiden_clus',
group_by = 'selections',
color_as_factor = TRUE,
point_size = 1,
coord_fix_ratio = 1,
cow_n_col = 2,
show_legend = TRUE)
Finally, you can plot the regions previously selected.
plotPolygons(my_giotto_object,
x = my_spatPlot,
polygon_name = 'selections')